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ABSTRACT 


paicutations of transient thermal stress levels were 
made for composite coatings of hollow cylinders. Stress 
analysis was done during a thermal cycle of successive 
heating and cooling periods. Various coating configurations 
were investigated. The results showed the effect of changes 
in coating arrangements and the failure modes due to the 
excessive stress levels. 

Solutions were obtained by approximate techniques, and 
a computer program for transient three dimensional thermal 
stress analysis was developed. The finite element technique 
was employed for the determination of the temperature 
distribution and the elastic stress analysis. Fortran IV 


G Level was used as the programming language. 
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PART I. THERMAL STRESS ANALYSIS OF 
CYLINDERS WITH COMPOSITE COATINGS 


The efficiency of a gas turbine engine can be increased 
by raising the turbine inlet temperature. However, protec- 
tion against high temperature oxidation or sulfidation of 
the turbine blades must be provided. One method, which has 
been suggested, [1]*, is the application of ceramic coatings 
which resist high temperature corrosion. Since the material 
properties vary through such a composite blade, thermal 
stresses become an important consideration particularly 
during transients with non-uniform temperature distributions. 

The application of these ceramic coatings has been 
accomplished via a process of high-rate sputtering [1]. 

This process requires high temperatures to insure proper 
adhesion of the coatings. Since the thermal life cycle of 
the coatings differs from the thermal conditions when the 
coatings are applied, the deposition temperature must be 
considered in the design criteria. 

In order to determine the thermal stress levels within 
possible coating arrangements, hollow cylinders were analyzed 
with various material compositions. A thermal cycle con- 
Sisting of a rapid heating phase followed by a cooling phase 


was used to generate temperature fluctuations and the thermal 


* Numbers in brackets refer to similarly numbered 
references in the bibliography. 
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stress distributions. The finite element technique was 
used to perform the elastic stress analysis and to determine 


the temperature distributions. 
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I. SIMULATION OF THE EXPERIMENTAL MODEL 


An Ep pecimentat program [1] to develop the coating 
technique and to observe the coating during thermal transients 
was begun using thin walled cylinders, a simple geometry 
instead of the more complex turbine blade. Such cylinders 
were prepared by the high-rate sputter deposition of nickel 
and stabilized zirconia on a nickel superalloy substratum. 

In addition to the pure coatings, a graded composition layer, 
composed of nickel and zirconia, was used between the sub- 
stratum and the outer layer. The cylinders were exposed to 
a thermal cycle consisting of a rapid heating phase followed 
by a cooling phase. The process was accomplished by the 
cross-flow of high temperature air past the hollow cylinder, 
followed by the use of low temperature air. The analysis 
reported herein attempts to analytically simulate this 
experiment. In doing so, several simplifying assumptions 
were made. 

The gas flow around the cylinder was considered symmetric, 
therefore, only one half of the hollow cylinder was modelled. 
A plane, along the cylinder axis, passing through the forward 
stagnation point and the after stagnation point was considered 
eneddtabatic surface. NO axial variation of the heat trans- 
fer conditions was allowed, however, the local heat transfer 


conditions were varied in the circumferential directions. 
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Although during an actual experiment the cylinder was 
supported by some apparatus, the modelled cylinder was 
unsupperted. Therefore, all transfer of heat occurred only 
in the hollow cylinder, and it was free to expand or contract 
in all directions. 

The simulation analysis was used to determine the tem- 
perature distribution and the stress distribution present 
in the cylinder at various times throughout the thermal 
cycle. Excessive stress levels were then used to predict 
when and what type of failure would occur. 

The classical thermal stress equations are coupled to 
the heat equation, however, during this simulation analysis, 
the heat equation was solved separately. The temperature 
distribution was then used to solve the thermal stress 


problem. 
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II. CYLINDER CONFIGURATIONS 


Four configurations were analyzed to determine the 
stress levels at various times in the thermal cycle. All 
geometries consisted of a hollow cylinder three inches long, 
an inside diameter of 0.5 inches, and an outside diameter 
dependent on the particular coating used. The thickness of 
the nickel superalloy substratum was 0.05 inches, and the 
thickness of the coatings varied from 0.005 to 0.009 inches. 


Drawings of all configurations are shown in Figure l. 


A. CONFIGURATION NUMBER ONE 

The first coating investigated was a 0.005 inch layer of 
zirconia on the nickel superalloy. Material properties 
changed from those of the nickel superalloy to those of the 
Ziecconia coating at 0.3 inches radius. 

The finite element mesh used in this arrangement con- 
sisted of 21 isoparametric quadratic elements totaling 200 
nodes. The element discretization was specified as one ele- 
ment in the axial direction, three elements in the circum- 
ferential direction, and seven elements in the radial direc- 
tion. The seven radial elements were proportioned into four 
elements through the substratum and three elements through 


the layer. 


B. CONFIGURATION NUMBER TWO 
The second arrangement analyzed consisted of 0.002 inches 


of pure nickel on the substratum followed by 0.005 inches of 
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stabilized zirconia. The material properties abruptly 
changed at two radial positions. 

The finite element mesh for this configuration consisted 
of 21 isoparametric quadratic elements using 200 nodes. The 
subdivision in the axial and circumferential directions 
remained as in configuration number one. The seven radial 
elements were arranged with three elements in the substratum 


and two elements in each layer. 


C. CONFIGURATION NUMBER THREE 

The third arrangement analyzed consisted of 0.004 inches 
of graded layer on the substratum. This graded layer allowed 
a gradual change of the material properties from those of 
nickel to those of zirconia. This arrangement was used to 
match the material properties of the substratum to those of 
the coatings. A 0.003 inch layer of stabilized zirconia was 
placed on top of the graded layer. Properties of the graded 
composition varied in a linear fashion as described in Section 
eT. 

This configuration was discretized into 24 finite elements 
totaling 226 nodes. The axial and circumferential directions 
were divided into one element and three elements respectively. 
The eight elements in the radial direction were arranged with 
four elements in the substratum, two elements in the graded 


layer, and two elements in the outer layer of zirconia. 
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D. CONFIGURATION NUMBER FOUR 

The fourth configuration was composed of the substratum 
followed by a 0.002 inch layer of pure nickel. A 0.004 inch 
coating of graded composition was placed on the pure nickel. 
A 0.003 inch coating of pure zirconia was used as the outer 
layer. 

The finite element mesh for this configuration consisted 
of 27 isoparametric quadratic elements using 252 nodes. The 
axial and circumferential directions were subdivided as in 
the other configurations. Nine elements were used in the 
radial direction, of which three were placed in‘the substratum, 


followed by two in each of the three layers. 


E. COMMENTS ON THE FINITE ELEMENT MESH 

Emphasis has been placed on the radial and tangential 
variations, therefore, the mesh was made finer for those 
directions. Since a limitation on the computer core require- 
ment and computing time dictates the number of elements and 
nodes, the axial direction was modelled using only one element. 
The problem solution is a three dimensional method, there- 
fore, ideally it would be more appropriate to have additional 
elements. Since the main interest was focused on the radial 
variation in material composition and the circumferential 
variation of heat transfer coefficients, the mesh used was 
considered appropriate. There was no axial variation in 
material properties or boundary conditions, therefore, no 


axial variation in the temperature distribution resulted. 
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A slight variation in the stress distribution resulted due 


to the three dimensional nature of the problem. 
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(*) Radial dashed lines indicate element boundaries in the 
circumferential direction. 
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CYLINDER CONFIGURATIONS 
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Tit. MATEREAT ‘PROPERTIES 


The materials used in the analysis were a nickel-super- 
aeloy. (IN-738), pure nickel, and stabilized zirconia (Zr 02 
plus Ca 0). The substratum of the hollow cylinder was com- 
posed of the nickel superalloy, and the coatings were com- 
binations of nickel and zirconia. Material property values 
were assumed constant throughout the temperature range inves- 
tigated. The values used were the average of those available 
in the literature [2,3]. The numerical values for elastic 
modulus, Poisson's ratio, coefficient of linear expansion, 
thermal conductivity, and volumetric heat capacity (density 
times specific heat) of the various materials are given in 
Table l. 

Matching of material properties was attempted by using 
a graded composition between layers in some of the coating 
configurations. Approximation to the properties of these 
compositions was made by Higcing a function based on the 
radial position. If two materials A and B have properties 
P. and Pp the property of the composition P, was determined 


A G 
by 


Se) (1) 


where rn is the radius of material A and Ine that of material 


B. 
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TABLE 1 


MATERIAL PROPERTIES 


Maekels Superalloy = IN-738 


Modulus of Elasticity 
Poisson's Ratio 


Coefficient of 
Linear Expansion 


Thermal Conductivity 


Volumetric Heat Capacity 


Pure Nickel 
Modulus of Elasticity 
Poisson's Ratio 


Coefficient of 
Linear Expansion 


Thermal Conductivity 


Volumetric Heat Capacity 


Stabilized Zirconia 
Modulus of Elasticity 
Poisson's Ratio 


Coefficient of 
Linear Expansion 


Thermal Conductivity 


Volumetric Heat Capacity 


31.75 x 10° PST 


0.3 


8.6 x Tom” in/in-deg.F 


2 


1.34 x 10 “© BTU/min-in-deg.F 


02029 BIU/cu.in.—deg.F 


30.00 x 10° psi 


0.3 


SG x ome in/in-deg.F 


2 


6.89 x 10 © BTU/min-in-deg.F 


02035 BIU/cusin.-—deg.F 


50h x om in/in-deg.F 


3 


3.36 x 10 ~ BTU/min-in-deg.F 


0.0382, BTU/ cu. in.=-deg.F 
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The materials used in the graded layer were pure nickel 
and stabilized zirconia. The values for the various pro- 
perties of the graded compositions are shown in the material 
property graphs, Figure 2. The property values used during 
the quadratic variations were determined by specifying 
property values for three points in the layer and then forming 
the quadratic function. Two of the points were specified 
to have properties of the inside and outside of the graded 
layer, i.e., the properties of nickel and zirconia. The 
midpoint was specified as the arithmetic mean of the pro- 
perties at the midpoint of the linear variation and the 
properties of zirconia. The resultant quadratic variation 


is shown on the same plots as the linear variation. 
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FIGURE 2. 


MATERIAL PROPERTY GRAPHS 
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IV. THERMAL CYCLE CONDITIONS 


The conditions which specify the thermal stress problem 
are classified as boundary conditions and initial conditions. 
With respect to the temperature distribution problem, the 
boundary conditions consist of convection heat transfer 
coefficients and the external fluid temperatures, and the 
gutecital condition is the initial temperature distribution. 
The stress problem has boundary conditions which specify 
zero or given displacement at various positions. In addi- 
tion, the zero-stress temperature distribution must be 
supplied since all temperature distributions are compared 


to the zero stress condition. 


A. HEAT TRANSFER BOUNDARY CONDITIONS 

The thermal cycle consisted of a rapid heating phase and 
a cooling phase. The cylinder wall was heated from room 
temperature (77 deg. F) to approximately 2000 deg. F ina 
period of 45 seconds. During the cooling phase, the heat 
transfer conditions were changed to cool the cylinder back 
to room temperature in approximately two minutes. During 
the heating phase, a convection boundary condition was main- 
tained on the external surface and an adiabatic condition on 
the internal surface. The cooling phase was characterized 
by convection boundary conditions on both internal and 


external surfaces. A circumferential variation of the 
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external heat transfer coefficients was used during both 
phases. The internal heat transfer coefficient was uniform 
during the cooling phase. No axial variation in the heat 
transfer coefficients was used. 

In order to determine an average heat transfer coeffi- 
cient for both the heating and cooling phases of the transi- 
ent cycle, the cylinder was treated thermally as a lumped 
capacitor. Using this average coefficient, an average 
Nusselt number was calculated using air as the external 
heating fluid. The temperature of the external fluid was 
assumed to be 2500 °F. A Reynolds number for flow past a 
cylinder was then computed using the Hilpert correlation. 
This Reynolds number was used to determine the local Nusselt 
number and therefore heat transfer coefficient for air flow 
past a horizontal cylinder. These values of Nusselt number 
were obtained from Kreith [5] and Krall and Eckert [6]. 

For the cooling phase, the average heat transfer coefficient 
was used for the internal surface. The values used during 
the computer analysis are shown in Table 2. These values 
were used for all configurations. 

The resultant temperatures of the cylinder wall were 
considered adequate for the analysis. The temperature 
variation at 45 seconds was between 1850 deg.F and 2050 deg F. 
Room temperature (77 deg. F) was reached after 2.25 minutes 


into the cooling phase. 
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TABLE 2 


Heating Phase 


Piura Temperature: 2500 deg. F 3 
Average Heat Transfer Coefficient: 22.6 BTU/hr-ft*~-deg. F 
Average Nusselt Number: 22.5 

Reynolds Number: 2262 

Equivalent Fluid Velocity: 124.6 ft/sec 


Local Heat Transfer Coefficients 


Angle Position from Stagnation Coefficient 
(degrees) (BTU/hr-ft2-deg. F) 
0.0 40.0 
30.0 41.0 
60.0 38.0 
510.0 26.0 
12:0'0 pois 5 
150.0 16.0 
180.0 i250 


Cooling Phase 


Fluid Temperature: 60 deg. F 3 
Average Heat Transfer Coefficient: 13.5 BTU/hr-ft' -deg. F 
Average Nusselt Number: 47.1 

Reynolds Number: 7300 

Poui1valent Fluid Velocity: 22.9 ft/sec 


Local Heat Transfer Coefficients 


Angle Position from Stagnation Coefficient 
(degrees) (BTU/hr-ft2-deg. F) 

0.0 22.0 

30:0 20.25 

60.0 15-0 

90-0 95.0 
22°0:.:0 6.0 
P5020 27.0 
180.0 ¥220 
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B. STRESS ANALYSIS BOUNDARY CONDITIONS 
The hollow cylinders analyzed were free to expand or 


Gontract in all directions. 


C. INITIAL CONDITIONS 

The cylinder was specified to be at an initial uniform 
temperature of 77 deg. F. The zero-stress temperature 
distribution was assumed to be the temperature when the 


coatings were applied. This was specified to be 932.0 deg. 
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VV; RESULTS..OF THE. ANALYSIS 


The results of the analysis showed that the stresses 
within the surface coatings varied from compression at room 
temperature to tension at the maximum temperature. Due to 
the variation of the heat transfer boundary conditions, a 
uniform temperature equal to the zero-stress condition was 
not present during the thermal cycle. The complex nature 
of the circumferential variation in heat transfer coeffi- 
cients causes temperature and stress distributions vastly 
different from those common to radial heat flow. 

All configurations follow the same general stress history. 
During the period when the temperatures are below the zero- 
stress condition, the coatings remain in compression. At 
the same time, the sign of the substratum stresses varies 
according to the position. At approximately 12 seconds, the 
temperature distributions reached the zero-stress condition 
and the stresses in the coatings became tensile. These ten- 
Sile stresses increased as the temperature increased 
reaching a maximum at 0.75 minutes. All configurations 
followed this general scheme. The maximum tensile stress 
is dependent on the zero-stress condition, i.e., the higher 
the zero-stress temperature, the lower the maximum tensile 
stress. Configuration Four had somewhat lower stress levels 
than the other configurations; however, the difference is 


not considered significant. 
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The use of a graded layer was investigated further by 
specifying a quadratic composition instead of a linear 
composition. An adjustment of the thermal conductivity and 
the coefficient of linear expansion to follow a quadratic 
function produced only minor changes in the stress levels. 
When all properties were adjusted to the quadratic composi- 
tion variation, no appreciable change occurred in the 
stress levels. 

Although the maximum stress levels are excessive, the 
maximum strains (9.6 x 1:0; >) are considered small enough to 
justify the use of a linear elastic analysis. In addition, 
the brittle nature of zirconia justifies the use of a brittle 
fracture failure criteria. Therefore, the elastic analysis 
provides a reliable method for determining whether the struc- 
ture will fail or not. Using the maximum principal stress 
theory of failure, some prediction of the time and tempera- 
ture when failure occurs can be made. If the ultimate 
strength of the ceramic material is set at 18.0 KSI and 
that of the nickel superalloy at 90.0 KSI, observation of 
the average tangential and axial stress levels for all con- 
Exguracivons indicates that failure would occur at about 0.2 
minutes and an average temperature of 1000.0 °F. At most, 
enas is only a rough calculation. Observations of an actual 
specimen could confirm this estimate by raising and lowering 
the specimen to successively higher temperatures, followed 


by appropriate microscopic examination. 
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The graphical representations of the thermal stress 
and temperature distributions were made using the numerical 
values obtained for the nodes at the cylinder ends. No 
attempt was made to smooth out the curves and no interpola- 
tions were performed. Each plot is identified by a time and 
a configuration designation. For example, the first plot is 
identified by TIME = 0.05 CONF.1-0. This plot represents 
the conditions at 0.05 minutes after start of the thermal 
cycle for configuration number one at position zero degrees 
relative to the stagnation point. Each axis has been labelled 
and the units identified. Since the range of values for 
stress and temperature varies from one graph to another, the 
scale factors used in plotting are provided on each plot. 
In order to provide an overall view of the stress levels 
throughout the thermal cycle, a plot of extreme principal 
stresses versus time has been prepared for each configuration, 


eraures 7, 14, 21, 26. 


A. CONFIGURATION NUMBER ONE 

This configuration contained a single layer of stabilized 
Zirconia on the nickel superalloy substratum. Comparison 
of the maximum principal stresses in all four configurations 
shows that configuration number one has the highest stress 
level. This maximum stress level (140.0 KSI) occurred at 
the maximum temperature, i.e., at the end of the heating 


phase, and was located at a radius of 0.3017 inches, the 
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zero angle position, and at the midpoint of the cylinder 
length. The largest tensile stresses present in configura- 
tion number one were all tangential and axial components 

of the stress state. The following figures show the radial 
variation of temperature and stresses at a given angular 
position for selected times during the heat transfer cycle: 
Figure 3 (a through j) at zero degrees, Figure 4 (a through j) 
at 60 degrees, Figure 5 (a through j}) at 120 degrees, Figure 6 
(a through j) at 180 degrees. (Similar sets of figures are 


provided for the remaining three configurations.) 
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B. CONFIGURATION NUMBER TWO 

This configuration contained a layer of pure nickel 
between the substratum and the layer of zirconia. The 
Maximum principal stress (134.9 KSI) occurred at the maximum 
temperature condition and was located at a radius of 0.3045 
inches, the zero angle position, and at the midpoint of the 
cylinder length. As with the other configurations, the 
larger tensile stresses were in the tangential and axial 
directions. 

The comparison of stress and temperature distribution 
as a function of angle position is shown in Figure 12. The 
axial stress level is nearly constant as the circumferential 
position changes, however, the radial and circumferential 
stresses change at about 90 - 120 degrees. This change in 
the stress level coincides with a change in the temperature 
distribution. The change in the temperature distribution 
results from the changes in the temperature boundary conditions 
in the down stream portion of the air flow past the cylinder. 
Figure 12 was prepared using data for the cylinder ends at a 
Badius of 0.3045 inches and a time of 0.75 minutes. 

The time dependence of the stress levels is shown in 
Figure 13. The tangential and axial stresses for a position 
on the cylinder end at a radius of 0.3045 inches and an angle 
of 60° vary in the same manner as the temperature at that 
point. The radial stress is minimal throughout the time 
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C. CONFIGURATION NUMBER THREE 

This arrangement of coatings included a graded layer of 
nickel and zirconia. The effect of this layer was to 
eliminate the abrupt change in the stress levels seen in the 
first two coating arrangements. The maximum principal stress 
(134.9 KSI) occurred at the maximum temperature. The minimum 
Peamcipal stress (-97.2 KSI) occurred at the beginning and 
end of the thermal cycle, i.e., the room temperature condition. 
Since the tangential and axial stress levels are excessive 
for a Major part of the cycle, failure would occur due to 
these stresses. 

As with configuration two, a comparison of stress and 
temperature distributions as a function of angular position 
was prepared using data for the cylinder ends at a radius 
@260-3055 inches and at a time of 0.75 minutes. The graphical 
representation of this comparison is shown in Figure 19. The 
function of stress and temperature versus time for the 60° 
angular position, at a radius of 0.3005 inches at the cylinder 


ends is shown in Figure 20. 
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D. CONFIGURATION NUMBER FOUR 

This coating arrangement included a layer of pure nickel 
between the substratum and the graded composition of 
configuration number three. The maximum principal stress 
(130.3 KSI) which occurred at the maximum temperature 
condition, was located at a radius of 0.3075 inches, the 
zero angle position, and at the midpoint of the cylinder 
length. As with the other configurations, the tangential 


and axial stresses were the largest for this configuration. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The thermal stress analysis of the various cylinder 
configurations showed that the maximum stress levels of the 
four coating arrangements were approximately equal. The 
maximum stress level was excessive enough to cause failure 
of the coatings and was dependent on the deposition temp- 
erature and the maximum temperature distribution. The abrupt 
change in the stress levels between the substratum and the 
outer coatings of configurations one and two was reduced by 
the use of the graded layers in configurations three and 
four. 

The following recommendations are suggested for future 


development and analysis. 


A. MATERIAL PROPERTIES 

An extensive investigation of the temperature dependent 
properties should be completed. In addition to the various 
individual materials, the graded layer properties must be 
thoroughly understood before a complete design analysis 
can be completed. All creep data for these materials must 
be obtained. If possible, material properties for the 


compressive stress state should be obtained. 
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B. VARIATIONS TO THE PROBLEM 

The following situations should be investigated. 

1. Other geometries, such as wedges and turbine blades, 
would exhibit different temperature and stress distributions. 
Investigation of these geometries may show lower stress 
levels which would be less susceptible to failure. These 
investigations should be made. 

2. The temperature distribution and maximum stress levels 
are directly related to the thickness of the zirconia coating. 
This coating depth should be varied to determine an optimal 
thickness for the minimization of the stress levels. 

3. Since the deposition temperature affects the maximum 
tensile stresses, more investigation to determine the most 
appropriate temperature for coating application should be 
made. More information about the compressive state will 
also allow a better evaluation of the effect of the deposition 
temperature selection. 

4. Although changes in the material property variation 
of the graded layer produced minor changes in the stress 
levels, the effect of the thickness of the graded layer should 
be investigated. 

5. Since. thermal load.conditions. are only part of the 
total load conditions on a turbine blade, an investigation 
of the effect of body forces on the stress levels in the blade 
should be completed. The program used in this analysis allows 


the specification of these transient body forces. 
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PART II. TRANSIENT THERMAL STRESS ANALYSIS 


The classical equations governing thermal stress analy- 
sis are coupled to the heat conduction equations, however, 
in most studies of this topic including this one, the stress 
analysis and the heat conduction are separated. The analysis 
of the transient thermal stress condition is performed in 
two steps. In the first one, the unsteady heat equation is 
solved for the temperature distribution. Then, using this 
temperature distribution, the equations of elasticity can be 
solved to determine the stress distribution. Analytically 
this solution process is possible only in the simplest of 
problems. However, through the use of the finite element 
techniques, an approximate solution can be obtained for more 
realistic problems. 

The computer program which was used as a basis to solve 
the heat conduction equation was written by Lew [7]. Several 
modifications were made to this program as described in 
Section I. The computer program used as a basis to solve 
the stress problem was written by Leonidas [8]. Several 
changes, as described in Section II, were also made to this 
program. 

Since three dimensional finite element analysis generates 
such a voluminous amount of information, an additional post- 
processing program was written to condense the relevant 


information. 
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All programs were written using Fortran IV G Level 
language and solutions were obtained using the IBM 360/67 


computer system. 
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I. TEMPERATURE DISTRIBUTION 


The unsteady heat conduction equation, expressed in 


rectangular coordinates, is 


om oT Boge ot zoe —) + = — 
me) * 5 5(Ky 55) tb aetk, gz) + 2 = pe (2) 


This equation is subject to the following boundary conditions: 


r= T, on the boundary or 
oT oT ow = 
hse ax 1, te i Dy ly + ke na 1, ate 6 + h(T Te) = 0 


(3) 


on the boundary 


and a specified initial temperature distribution. The 


symbols used in these equations have the following definitions: 


Q - rate of internal energy generation per 
unit volume 


) - density 
C - specific heat 


PE gis = thermal Conductivity in the x,y, and Zz 
directions 


oe ee - the direction cosines at a point on the 
Y boundary 


- external heat flux 
- surface heat transfer coefficient 


Tr - fluid temperature 
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As described in Ref. 9, the unsteady heat equation can 
be transformed into a discretized form by the finite element 
technigue. Using the symbol [ ] for ann xn matrix, 


_ereror a row vector, and { } for a column vector, 

Erie elt} + fF} = 10 (4) 
where the temperature T of equation (2) has been expressed 
as the product of the nodal shape functions and the nodal 


temperature values. 


TiN 1) (5) 


The time derivative of the temperature =" has been expressed 
as 
Pe x ° is e 
ae P= <N>) {Tt (6) 
The matrices [H], [C], {F} are defined at the element 
level by 


¥ ae oN 
a. = heer (a ~}< i= Sie a Ko )dx dy dz 
= (7) 
[Cl]. = f{N} pC <N> ax dy dz (8) 
V 
e 
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{F}, = - SQ{N}dx dy dz + sq{N}dS, + ( f {N}h<N>dS_) {T, } 
V S 
e S (9) 


where 
QF Bey h T, 


The program developed by Lew [7] solves the discretized 
equation subject to the listed boundary conditions. In order 
to adapt the program to the thermal stress problem, several 


alterations were incorporated. 


A. TIME INTEGRATION SCHEME 

Reference 7 provides numerous time integration schemes 
to deal with the first time derivative of the temperature. 
In this study only the best integration schemes were retained. 
A five point backward difference formula was used to approxi- 
mate the time derivative. For the ith time step, this 


approximation can be expressed as 


fo fo. = (3r ale i  Oh 


24 3 cre 


= PALS iH kane) 
al 


2 1 


(10) 


where At is the specified time step. The truncation error 
: ee he eta tT 
for this expression is +=(At) anor 
z ot 
If equation (10) is written in vector form and substituted 
into equation (4), the final form of the discretized heat 


equation is obtained. 
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me2at(H) + 25[c])-{T, } = 


peree,? ~ sett. 4) + W6eiT, ,} - 3{T,_,} - l2aetr} 

Ci) 
meni), (C], {fF}, At and the four preceding values of {T, } 
are known, the system of simultaneous equations is solved 
mor the next set of {T, }. 

The initial four values for the temperature distribution 
is obtained by using the trapezoid method for a specified 
number of steps at a small At. After each n steps at At, 
the values of tT, } are saved. When three sets of ti} plus 
the initial known distribution are available, the transient 
solution using a time step of (nAt) can proceed. For exam- 
ple, an initial step size of 0.001 taken for 150 steps 
results in a step size of 0.05 for the remainder of the 
transient problem. 

Using the solution technique described above, the trans- 
ient temperature distribution for an infinite hollow cylinder 


were determined. The problem characteristics are: 


Inside diameter = 0.5 inches 

Outside diameter = 0.6 inches 

Initial temperature =~ 932 Geq.WF 

Internal and external 

heat transfer coefficients = 25.0) BLU nr—-SG.rt. Ged. 
Internal Fluid temperature =. 896 deg. °F 

External Fluid temperature —— LOO 4 deg... 

Initial time step size = 0.001 min. 


150 


Number of initial steps 
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Results of the approximate solution and the exact solution 


are given in Table 3. 


B. NODAL MATERIAL PROPERTIES 

The usual method of specifying properties to be constant 
throughout a finite element becomes troublesome when com- 
posite materials have continuously varying properties. Many 
separate elements would be required in order to adequately 
represent a structure which has material properties that are 
continuous functions of spatial coordinates. 

In order to alleviate the problem, several routines in 
the basic program of Ref. 7 were altered to accommodate 
nodal properties in addition to element properties. A 
property P is expressed as P = N; Pie where P; is a vector 
of nodal values and Ns are the same shape functions used in 
equation (5) for the temperature variable. The matrices 
[H], [C], and {F} are formed by the same integration tech- 
niques as used in the basic program. However, the proper- 
ties are no longer constant during the integration, and 
must be re-evaluated at each integration point. 

As an example, consider steady state heat flow through 
wmeylinder in which the thermal conductivity varies contin- 
uously along a radius. As proposed above, the development 
of a twenty nodal point element could represent exactly any 
linear or quadratic variation of material properties. Table 
4 displays the results of such a problem as compared to 


the exact solution. The column heading "Element Property 
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TABLE 3 


TRANSIENT TEMPERATURE OF THE INFINITE HOLLOW CYLINDER 


Exact Solution (Deg. F) 


Radius = — 0.25 O25 10) 3.0 

(in) (min) 

0.25 12775165 1409.796 1479 .776 1491715 
0.26 1278 .061 1410.916 1481.013 VAG2 997 
0.27 1279..104 1412 .062 1482 .214 1494.182 
0.28 1280;,291 1413.236 1483.381 1495.348 
0.29 1281.618 1414, 438 1484.518 1496.473 
0.30 1283.081 1415.669 1485.626 1497.561 

Finite Element Solution 
Using 5-Point Finite Difference Method (Deg. F) 

Radius t= 0625 0.5 0 3.0 

(in) (min) 

0.25 12772220 1409 .897 1479. 879 1491.810 
0.26 IZ 76. VG 1411.016 1481.116 1493.068 
On27 1279 .160 1412 .163 1482 .316 1494.278 
0.28 1280 .347 1413.336 1483.483 1495. 443 
0.29 1281.674 1414.538 1484.612 1496.568 
0.30 293138 1415. 769 1485. 726 1497.654 


au 


4 
+ 
4 











Inside Diameter 


Outside Diameter 


Inside Surface Temperature 


Outside Surface Temperature 


TABLE 4 
STEADY HEAT FLOW THROUGH A HOLLOW CYLINDER 


0% 5="a'n. 


O.6 “an. 


[0050 °F 


1OOOSO °F 


Thermal Conductivity Variation in Outward Radial Direction: 
0505. to 0.25 BIU/in—-min-°F 


(a) 
(b) 


Radius 
Cin.) 


0.25 
Oe2625 
On275 
0.2875 
0.3 


Quadratic Variation of Thermal Conductivity 


Radius 
en. ) 


O25 
ORi2625 
0.275 
Oe 2875 
Or. 3 


EOr linear function, 


Fon quadratic function, 
(at mid-radius) to 0.25 BTU/in-min-°F 


Linear Variation of Thermal Conductivity 


(Two Radial Elements) 


Exact 
Solution 
ideg~ EF) 


500.0 
575.884 
669.145 
IDA. O97. 
1000.0 


Element Property 


Solution 


(deg. F) 


500.0 
5910,..807 1: 
677.910 
842.256 
1000.0 


(Two Radial Elements) 


EeaiGite 
Solution 
(deg. F) 


50050 
62'3..18 79 
748.213 
873.441 
1000.0 


Element Property 


Solution 


(deg. F) 


500.0 
627.602 
T49~253 
ST S42 
1000.0 


OBZ025 to 0.225625 


Nodal Property 
Solution 
(deg. F) 


500.0 
576.502 
67/07. 301 
MD6..19 8:2 
1000.0 


Nodal Property 
Solution 
(deg. F) 


5:0.0/.:0 
62:47..213 
748.803 
873: 747 
1000.0 
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Solution" identifies the results obtained by specifying 
properties constant over an element. The column heading 
"Nodal Property Solution" identifies the results obtained 
by specifying the properties variable over an element, 


in other words nodal properties. 


C. BOUNDARY CONDITIONS 

The load vector represented by equation (9) contains 
two surface integrals which pertain to boundary conditions 
described by equation (3). The development of this load 
vector is shown in Ref. 9. 

The usual technique for specification of ene noandeny 
conditions is to consider q and h as constants over the 
applicable element surface. In order to specify convection 
boundary conditions which vary over an element surface, the 
flux and heat transfer coefficient have been treated as 
nodal values in the same manner as nodal material properties. 
The heat generation term, Q, is also treated in this manner. 
Therefore, instead of the constant terms in equation (9), the 
usual shape functions N. are used with the nodal values to 
form the variable terms Q = N; Qi q = N, qa and h = N; h;. 
The integration processes proceed as described in Ref. 7, 
with the addition of these terms. 

The surface integrals are computed using a coordinate 
transformation to allow the use of Gauss Quadrature. This 


transformation maps the Cartesian coordinates (x,y,z) into 
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te 
v 


a local coordinate system (&,n,t). As described by Hanson 
[10], the unit area vector in the (x,y,z) system is evaluated 


in the (&,n,vc) system by 


a Ay i Seagrs 
n ds, = oem 3p) da débB , (1:25 
where 
me—ex(h.n, ct + ¥(E;net)3 + 2(E.n,t)k (13) 


Specifies the surface and a, represent —,n, Or t. 
The vector V = = x = is normal to the surface area 
da dg and must be evaluated during the integration process. 


It is convenient to express equation (12) as 
n ds. = |¥| ee da dg = [Vin do dg , (14) 


nw 
Mere Nn 1S a unit normal vector to surface area do dg. 


Therefore, during the integration process, the factor 


iw) = (2k 2% - 32 ayy” 4 (ae Bz _ Bz amy” 
da Of a dp da df da 0B 
De /2 
ax ay _ ay ax 
> Ns Sn ae Oe se 


must be evaluated. If the surface integral is performed 


Over an area on which one of the variables x, y, or z is 
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constant, the quantity |V| becomes identical to the deter- 
minant of the Jacobian matrix for the surface transformation. 
For a general surface r(x,y,z), the complete expressions for 


|v| must be used. 
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ER. | «SERESS DISTRIBUTION 


The finite element technique provides the required 
methods for the solution of the three dimensional elastic 
stress analysis problem. Reference 9 describes the deri- 
vations for the following equations. A structure can be 
divided into several elements, each formed by a specified 
number of nodes. If this structure is acted upon by various 
forces, the nodal displacements {r} are determined by 


solving the following equation: 


Pitz} = (k)} - (PF) - {Fh - tf), - (Fh, 6) 
O Oo 


where the matrices and vectors are described by 


[K] = teh) (D) [ey av (17) 
V 
ges ie ae 
as = ie [N] bes), av (18) 
Le Ju 
oe = - [N] {£3 av (19) 
fee -- -f, (s)" [p] fe.) av (20) 
oO V 
- At 
eo. = a PE to Jay: (21) 
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and {R} is the vector of concentrated nodal loads. The 


other symbols used in these equations are defined as: 


[N] = VeECtor Of shape functions 

[B] - matrix of shape function derivatives 
[D] - elasticity matrix 

{f,} - vector of body forces 

{£3 -- vector of surface forces 

fen) - vector of initial strains 

{o,} - vector of initial stresses 


Following the solution of equation(16), the nodal displace- 


ments {r} are used to evaluate the nodal strains. 


ten = EB) {xr} (22) 


The final step in the analysis is the calculation of the 


nodal stresses. 


cele lee = let) + tot (23) 


A computer program to make these calculations was 
written at the Naval Postgraduate School by Leonidas [8]. 
In that version of the program, the right hand side load 
vector was allowed to contain only concentrated loads {R}. 


In order to permit the presence of thermal stresses and body 
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forces, several changes were made to allow Cent and {£1} 


to have non-zero values. 


A. GENERATION OF LOAD VECTORS 

Temperature differences always generate strains but 
sometimes stresses may still be zero, for example, when 
the temperature differences are constant throughout, ¢€ in 
equation (23) would be equal to en giving zero stress for 
one material. The presence of a non-uniform temperature 
distribution in a solid body causes stresses and strains due 
to non-uniform expansion. In addition, a solid body composed 
of materials with differing coefficients of linear expansion, 
4, will at times experience stresses and strains. This happens 
when the temperature distributions differ from that of the 
zero-stress condition. 

The forces present in a body during a thermal stress 
state can be calculated by setting eie= AAT, where AT is 
the difference between the actual temperature distribution 
and the zero-stress temperature dilstri butksoneye Inhendere tec 
form a consistent load vector, these forces, calculated by 
equation (20), are added to any other nodal forces present 
on the right hand side of equation (16). 

The body forces, due to some force field, are converted 
into nodal forces by using equation (17). These are added 
Bovthe right hand vector of equation (16) to form the total 


consistent load vector. 
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B. NODAL MATERIAL PROPERTIES 

As discussed in Section IB, the use of nodal material 
properties is a more convenient method for composite body 
analysis. In the stress analysis operations, the elastic 
Matrix and the coefficient of linear expansion enter into 
the calculations. The elastic matrix is composed of the 
elastic modulus and Poisson's ratio. For the different 
operations, these properties are evaluated by using 
E = Nj Es, y= N, Use and A= N; Aye During the various 
integration processes these property values are evaluated 
at each integration step. 

The use of nodal properties requires careful considera- 
tion as to the effect the shape functions have on the pro- 
perties within the element. If an element of one material 
is adjacent to an element of another material, the use of 
nodal properties can lead to erroneous results. Since the 
nodes common to the elements can have only one set of pro- 
perty values, one of the elements will not have the expected 
uniform properties. It would be better to use constant ele- 
ment properties or combine the two elements into one and 
use variable properties over that element. The amount of 
property variation between the two elements will also assist 


the analyst in determining the proper technique. 
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fret. DEoCREPTION OF THE PROGRAMS 


Three programs were used to complete the thermal stress 
analysis. As noted in Sections I and II, the programs of 
Lew and Leonidas [7,8] were used as a basis to develop the 
programs for the temperature distribution and elastic stress 
analysis. The third program was written to process the 
information generated in the stress analysis. The computer 
listings for these three programs are provided in Appendices 


iG, and H. 


A. TEMPERATURE PROGRAM 

The main alterations to the temperature program [7] were: 
(1) the elimination of the second time derivative of the 
quasi-harmonic equation, (2) the addition of the time inte- 
gration scheme described in Section IA, and (3) changes to 
Subroutines FORMC, CUBE, and BCOND to allow the specification 
of nodal material properties and nodal boundary conditions. 
The version of the program listed in Appendix F allows 278 
nodes and a bandwidth of 40. These specifications were 
varied to adjust the computed core requirements to the 
particular problem being solved. A problem with 278 nodes 
and a bandwidth of 40 required approximately 350,000 bytes 
of storage. 

The assumptions made in the solution of a problem using 
this program are: (1) material properties are constant with 


time and independent of temperature, (2) directional material 
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properties must align with the global coordinate axes, and 
(3) the boundary conditions are constant within a heat 
transfer phase, but may be altered at the end. 
i Subroutine TIMEHB 

As discussed in Section IA, this subroutine performs 
a time integration using a five point backward difference 
method. During each heat transfer phase, boundary conditions 
are constant. Following applications of these boundary con- 
ditions to the appropriate matrices and vectors, the trape- 
zoid integration method is used to generate three starting 
values. Along with the initial temperature distribution, 
these three starting values provide the necessary information 
to use the five point difference scheme. The output of tem- 
perature distributions is controlled in a manner similar to 
that used in the original program [7]. These output values 
are adjusted to the reference temperature supplied in the 
data deck. During a stress analysis problem, this reference 
temperature is the zero-stress temperature. Since these 
temperature distributions are used in a subsequent stress 
analysis, proper control is provided to punch cards or 
eranster the information to temporary storage devices. The 
program in its present form allows up to twenty separate 
temperature vectors to be stored for future use. 

2. Subroutine FORMC, CUBE 
Subroutine FORMC was modified to allow nodal proper- 


ties to be used in the integration process for the matrices 
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{H] and [C]. The subroutine CUBE performs the cubic inte- 
gration used in the formation of these matrices and there- 
fore required slight modification. 
3. Subroutine BCOND 

Subroutine BCOND and its subsidiary routines were 
modified to accommodate nodal boundary conditions. Subrou- 
tine JACOB was altered to evaluate the factor given in 
equation (15). This factor was used in the surface integra- 
tion as discussed in Section IC. 

4. Data Preparation for the Temperature Program 

The data required for a transient solution is com- 
prised of eight groupings of cards. Careful detail must be 
followed to insure proper formats are used and the correct 
order of cards maintained. A data check can be performed 
without any major calculations. The data check normally 
requires 20-30 seconds, however, the full computer core 
requirement is needed, 

Preparation of the data cards is explained in the 
initial comment cards of the program listing in Appendix 
F. A sample data deck for a 32 node, 2 element program 
is provided in Appendix A. Several aspects of the data 
preparation must be thoroughly understood. 

a. Nodal Properties 

The method of input for values which are common 

to several nodes was provided to reduce the number of input 


Cards. Care must be taken to insure that every node has 
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properties assigned to it, and that the order of material 
property cards follow the order of ascending nodal numbers. 
b. Output Storage 
The number of times stored output is requested 
must match the number of loads used in the stress analysis 
program. 
ec. Steady State Solution 
If a steady state problem is being analyzed and 
no convection or internal generation conditions exist, a 
blank card must be inserted in place of the heat transfer 
phase data. 
ad Stacking of Problems 
Tr this program is used in conjunction with the 
elastic stress analysis program, only one problem can be 
loaded at a time. Therefore, the final comments concerning 
stacking problems do not apply, and a "STOP" card must be 
used. 
e. Boundary Condition Data 
When preparing the boundary condition data, the 
number of cards with boundary property values (NODECD) is 
given on the first card of each boundary condition data 
set. Ensure that NODECD cards are used in that particular 
set. 
5. Job Control Language 
The necessary job control language for using the 
temperature program is listed in Appendix B. The two 


@ards beginning with //GO.FT1OFO01 and //GO.FT11F001 are 
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mecesSsary for the storage of the matrices [H] and [C]) during 
the time integration process. They are not required for 

the steady state solution. The card beginning with 
//GO.FTO7FO01 is used to store temperature distributions 

on temporary files for use in the stress analysis program. 
If punched output is desired, this card should be replaced 
by f/GO.ETOJEOOL..DD: .SYSOUT=B 

The card beginning with //GO.FTO6F001 is necessary to allow 
large amounts of output from the computer printing devices. 
If less than 3000 lines of output is anticipated, this card 
can be removed. These job control cards are designed for 
use at the W.R. CHURCH COMPUTER CENTER, NAVAL POSTGRADUATE 
SCHOOL, Monterey, California. Use of the program at another 


installation will require appropriate modifications. 


B. ELASTIC STRESS ANALYSIS PROGRAM 

Several changes were required to add a transient thermal 
stress analysis capability to the stress analysis program 
[8]. These included: (1) Addition of a consistent load 
generation subroutine LODGEN, (2) Alteration of the subrou- 
tine SOLVE to allow multiple loads, (3) Changes to subrou- 
tines INPUT, SORT, STRESS to provide steps for various 
calculations pertaining to the temperature distributions, 
(4) The capability to use nodal properties requiring the 
addition of subroutine NDPROP and numerous alterations to 
the integration process for the computation of the "stiffness 
Matrix", and (5) The addition of a subroutine PRNSTR to 


Calculate the magnitude of the principal stresses. 
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The program uses external storage for all information. 
Therefore, numerous additions and changes were made to the 
various WRITE and READ statements. A thorough knowledge 
of the direct-access statements used in Fortran IV is 
essential to understand the program used in this analysis. 

1. Subroutine LODGEN 

This subroutine performs calculations used to form 
the consistent load vectors as described in Section IIA. 

\ 

A lengthy computing time is required for this subroutine due 
to the number of computations and the many direct access 
file operations. The number of Gauss points (NGPL) used 
can be specified independently of the number used for the 
formation of the "stiffness matrix". The consistent load 
vector forms the right hand side of equation (16) and must 
be of sufficient accuracy to insure valid final results. 
For the calculations made in Part I, four Gauss points were 
considered sufficient. 

2. Subroutine SOLVE 

The basic subroutine used to solve equation (15) 
is contained in Refs. 8 and 11. In order to solve multiple 
load vectors, changes were incorporated to make necessary 
computations to all load vectors simultaneously. As des- 
Cribed in Refs. 8 and 11, calculations are performed on 
each block of values in the load vector until all nodes 
have been processed. The changes which were made performed 
these calculations to a specified block for all load vectors 


prior to advancing to the next block. This method was used 
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to minimize the number of computations and the direct 
access operations. 
3. Subroutine NDPROP 

This subroutine calculates the material properties 
at a node or Gauss point. If the elastic matrix is re- 
quested, the various properties are calculated for the 
position in question, and the matrix is formed. If the 
coefficient of thermal expansion is requested, it is calcu- 
lated for the position specified and returned to the calling 
program. This subroutine uses a separate subroutine, SHAPE, 
to evaluate the shape function values used in the calculation 
of the property values. 

ie Preparation, of Input. Data 

The preparation of all data cards is explained in 
the initial comment cards of the program listing, Appendix 
G. Special care must be given to the details of the input 
cards. 

a. Object Time Formats 

The use of object time formats is made throughout 

the input subroutine. The use of object time formats elimin- 
ates the need to recall what format a certain data card is 
tO follow. The object time format technique is covered in 
most books on the Fortran language. 

b..4 Data check 

A data check can be made by specifying NGPS=0. 

thas requires the full 250,000 bytes of storage, but is 


completed in less than 20 seconds. 
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cc. External Storage of Load Conditions 

If the load conditions of temperature distribu- 
tion and body forces are prepared by a preprocessing program, 
these values are passed using files 33 and 34. File 35 for 
use with nodal displacements was added to the input routine 
for future development. When the values are passed from a 
preprocessing program, a time value must be provided at the 
beginning of each set of load conditions. 

5. Job Control Language 

Since the stress analysis program stores all infor- 
mation on external devices, the job control language is a 
very critical part of the program. A very thorough under- 
standing of the specific requirements for a particular com- 
puter facility is needed to use this program. All examples 
and instructions presented here apply to the W.R. CHURCH 
COMPUTER CENTER AT THE NAVAL POSTGRADUATE SCHOOL, Monterey, 
California. 

External storage operations used in the program are 
divided into two types. One type is the sequential file 
READ or WRITE statement. This method of external storage 
performs the operation on a specified record length in a 
sequential manner. The second type of operation is the 
direct access method. In this case, a specified record 
length is written into or retrieved from a given location 
meenin the file. In the first type of operation, no 


changes to the program language are necessary. The second 
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type of operation requires changes to the DEFINE FILE 
statements listed in the control program. These statements 
will vary according to the size of the problem being 
analyzed. 

Each file used by either a sequential type operation 
or a direct access operation must be described in the job 
control language. For the case of the direct access opera- 
tions, the DEFINE FILE statements and the job control language 
must agree both in record length and the number of records. 

In order to provide a complete summary of all files, 
Table 5 was prepared. If a particular problem does not 
require a file, its job control language and any DEFINE FILE 
statements can be eliminated. Table 6 describes each DEFINE 
FILE statement and job control card. For each file, a for- 
mula is supplied which determines the record length and the 
number of records which that file will contain. The program 
fescing, Appendix G, and the job control card listing, 
Appendix D, are examples for use with a thermal stress 


problem uSing 278 nodes, 27 elements, and 10 load vectors. 


©. POST PROCESSING PROGRAM 

The elastic stress analysis program provides a consider- 
able amount of information. For a problem of 200 nodes 
weeh 10 load conditions, about 16,000 lines of output are 
generated. Interpretation of this amount of output is both 


tedious and time consuming. Usually only a few nodes are 
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TABLE 5 


PLLE -DEEINETIONS 


File Number = Definition 
7 = Stotage of Total Stiffness Matrix Blocks 
8 - Storage of Element Stiffness Matrices 
9 - Storage of Element Connectivity 
10 - Storage of Nodal Coordinate Vectors 
ie - Storage of Nodal Properties 
EZ - Storage of Nodal Stress Vectors 
13 - Storage of Nodal Strain Vectors 
14 - Storage of Concentrated Load Vectors 
-5 - Storage of Boundary Condition Vectors 
16 - Storage of Element Material Numbers 
Li - Storage of Reaction Vectors 
18 - Storage of Nodal Temperatures 
ins) - Storage of Nodal Coordinates by Element 
20 - Storage of Load Vector Blocks 
7A - Storage of Nodal Displacement Vectors 
22 - Storage of Predefined Boundary Displacements 
23 - Storage of Nodal Properties by Element 
24 - Storage of Nodal Temperatures by Element 
25 - Storage of Nodal Body Forces 
26 - Storage of Nodal Displacements 
PAG) - Storage of Nodal Body Forces by Element 
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File Number 


28 


2) 
30 
33 


34 


35 


Definition 


Storage of Nodal Displacements 
by Element 


Storage of Nodal Force Vector 
Storage of Time Vectors 


Storage of Preprocessed Temperature 
Distributions 


Storage of Preprocessed Body Force 
Distributions 


Storage of Preprocessed Nodal 
Displacement Distributions 
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TABLE 6 


DEFINE FILE Statements and Job Control Cards 


General format for these cards: 


DEFINE FILE XX(NRX,LRWX,U,I1XxX) 
//GO.FTXXFO01 DD UNIT=SYSDA,SPACE= (LRBX, (NRX,1)) , DISP=(NEW, DELETE) 


where XX,NRX,LRBX are defined in the table. 


FILE Number of Record length Record length 
Records in BYTES in WORDS 
XX NRX LRBX LRWX 
1) MMxNNXNREC7 NCOUNTx8 /NREC7 NCOUNTx2/NREC7 
8 4xNEL 7200 1800 
9 NEL 80 20 
10 NDPT 24 6 
ab NDPT 24 6 
a2 NDPTxXNLOAD 56 14 
Ss NDPT 56 14 
14 NCLD 32 8 
5 NPBC 16 4 
16 NEL 4 a 
Ey NPBC 32 8 
18 NDPTXNLOAD 8 2 
19 NEL 480 120 
20 NNXNLOAD NSx8 NS x2 
2A NNxXNS/NDF 24 6 
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XX NRX LRBX LRWX 
22 NPBC2 24 6 
23 NEL 480 120 
24 NLOADXNEL 160 40 
25 NLOADXNDPT 24 6 
26 NLOADXNDPT 24 6 
Di NLOADXNEL 480 120 
28 NLOADXNEL 480 120 
29 NLOADXNDPT 24 6 
30 NLOAD 8 2 
33 NLOADXNDPT 8 2 
34 NLOADXNDPT 24 6 
35 NLOADxNDPT 24 6 
Notes: 
fay) For a 20 nodal point brick, NREC7 = 4 
(Bb) For a 32 nodal point brick, NREC7 = 6 
(c) Variables MM, NN, NS, etc. are defined in Ref. 8 
232 





considered at one time and the remainder of the information 
is reserved for later use. 

Through the use of job control language, the important 
information can be stored on external storage devices and 
processed by a separate program. A post-processing program 
was written to analyze the stored information. 

fe. Program Description 

The program listed in Appendix G performs the 


following operations: (a) information stored on external 
devices is read for a specified number of nodes; (b) the 


x,Y, and z components of all stress vectors are analyzed 
for the principal stress values, which are used to determine 
¢he maximum and minimum values for a given load condition; 
(c) Offline and Calcomp plots of stress versus coordinate 
position and temperature versus coordinate position can be 
obtained; (d) Offline and Calcomp plots can be obtained for 
the stress values in polar coordinates. 
2. Input Data Preparation 

The input data deck has been designed to consist of 
three or four cards. The preparation of these cards is 
explained in the initial comments of the program in Appendix 
H. When more than one problem is processed, care must be 
taken to limit the number of Calcomp plots requested. All 
Calcomp plots are performed twice to make the lines dark 


enough for reproduction by other means. 
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3. Job Control Language 


Since the information is to be passed from the 
stress analysis program, a few changes to the job control 
language discussed in Section IIIB are necessary. Files 
fopee2, 18, and 30 are changed to the following format: 


//GO.FTXXFO01 DD UNIT=XXXXX,SPACE=(LRBX, (NRX,1)), 
// DISP=(NEW,KEEP) , DSN=YYYYY, VOLUME=SER=Z2Z222Z 


where XX, LRBX, NRX - defined in Section IIIB 
XXXXX - a unit number 
WYYYY - a file name, e.g. STRESS 
ZZZ2Z - a device number 


The unit number and device number are obtained from the 
Barticular computer facility. 

The job control language for the post-processing 
program must match with that passed by the elastic stress 
Eialysis program. Files 8, 9, 10, and 11 are used for the 
stress values, coordinate values, temperature values, and 
time values respectively. The Jeb Control cards for the 
post-processor have the same form as those in the elastic 
stress analysis program. The linkage between the two programs 
is formed by specifying the same unit number, device number, 
and file name for the information given by the stress program 
and received by the post-processor. For example, file 12 
in the stress program would be 


@7GO.FTI2F001 DD UNIT=2311,SPACE=(LRBX, (NRX,1)), 
// DISP=(NEW,KEEP) , DSN=STRESS , VOLUME=SER=SYS001 


and, therefore, File 8 in the post-processor would be 
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//GO.FTO8FO01 DD UNIT=2311,SPACE=(LRBX, (NRX,1)), 
// DISP=(OLD,KEEP) , DSN=STRESS , VOLUME=SER=SYS001 ° 


The post processing program also contains DEFINE FILE 
statements which follow the guidelines in Section IIIB. 
Further details on the procedures for job control 


language can be found in Ref. 12. 


D. SYSTEM USAGE OF ALL PROGRAMS 
The discussions thus far have pertained to the individual 
programs, however, the actual production usage of the pro- 
grams was done in a multistep process. The usage of the 
programs as a system first requires the formation of the 
system, then the production run of a problem. 
1. System Formation 
Appendix D lists the job control language used to 
create the library of programs. Details of this procedure 
ere provided in Ref. 12. 
Zee eroduction Runs 
Appendix E lists the job control language for the 
multistep process of temperature and stress distribution 
determination. The major change in the methods is the 
transport of the temperature distribution from the first 
Step to the second step. This is accomplished by the card 
beginning with //GO.FTO7F001 for the first step job control 
language and the card beginning with //GO.FT33F001 of the 
second step job control language. These two cards connect 


the informatin generated and used by the two programs. 
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IV. RECOMMENDATIONS 


As with any endeavor to construct a usable tool for 
engineering analysis, the job is never done. The following 


changes are recommended for future development: 


A. PROGRAM FOR THE TEMPERATURE DISTRIBUTION 

The present version of the temperature distribution is 
designed for an "in-core" solution technique, which restricts 
the problem size to within the computer core size. The 
stress program uses an "out of core" technique and is 
restricted only by the availability of external storage 
devices. In order to eliminate the core size restriction, 
the temperature program requires complete revision to permit 


the use of external storage devices. 


B. ENETIAL STRESS DISTRIBUTION 
The load generator should be altered to allow a non-zero 


eieutial Stress distribution. 


C. JOB CONTROL LANGUAGE 

The amount of job control language allowed at a par- 
eveular Computer facility is limited. For the facility 
used during this analysis, it appears that no more job 
control language can be used with this program. The various 
files used in this program must be combined to reduce the 


humber of job control cards needed. 
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D. MATERIAL PROPERTIES 

In order to more accurately solve the real engineering 
problems, the capability of temperature dependent properties 
should be added to both the temperature program and the 


stress program. 
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APPENDIX A 
SAMPLE DATA DECK FOR TEMPERATURE PROGRAM 
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APPENDIX C 
SAMPLE DATA DECK FOR STRESS PROGRAM 
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APPENDIX D 
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